function [oem_profit, oem_profit_sim, carrier_ind_l, oem_ind_l, model_id_l, p, w, share,...
    cons_surplus, cons_surplus_no_error, expmu, carrier_profit, p_sim, ...
    w_sim, share_sim, prod_profit_binary, elas_ns, craft_l] = ...
    equi_ind(config, binary_prod, fixed_prod, sim_pcoeff, setup, ...
    para_nonlinear, randomdraw, ...
    rand_income, income_norm, expansion_nb)
% PURPOSE: compute pricing equilibrium and variable profit

n_oem = length(binary_prod{1}.owner_ind);
n_x_r = length(binary_prod{1}.x_rand);
n_x_r2 = length(binary_prod{1}.x_rand2);
n_prod_binary = sum(config);

if nargin<15 || isempty(expansion_nb)
    n_draws = size(binary_prod{1}.demand_resid, 2);
else
    n_draws = size(binary_prod{expansion_nb}.demand_resid, 2);
end

w0 = nan(n_prod_binary, 1);
p0 = nan(n_prod_binary, 1);

meanutility_l = nan(n_prod_binary, n_draws);
mc_mean_l = nan(n_prod_binary, n_draws);

oem_ind_l = nan(n_prod_binary, n_oem);
carrier_ind_l = true(n_prod_binary, 1);
model_id_l = nan(n_prod_binary, 1);
x_rand_l = nan(n_prod_binary, n_x_r);
x_rand2_l = nan(n_prod_binary, n_x_r2);
craft_l = false(n_prod_binary, 1);
ind_endo_l = true(n_prod_binary, 1);

k = 0;
for nb = 1 : size(config, 1)
    if config(nb)==1
        new_prod = binary_prod{nb};
        ind = k+1;
        
        w0(ind) = new_prod.wp;
        p0(ind) = new_prod.p0;
        
        oem_ind_l(ind, :) = new_prod.owner_ind;
        
        if ~(nargin<15 || isempty(expansion_nb)) && nb==expansion_nb
            meanutility_l(ind, :) = new_prod.demand_resid;
            mc_mean_l(ind, :) = new_prod.mc;
        else
            meanutility_l(ind, :) = new_prod.demand_resid + zeros(1, n_draws);
            mc_mean_l(ind, :) = new_prod.mc + zeros(1, n_draws);
        end
        
        model_id_l(ind,:) = new_prod.id;
        
        x_rand_l(ind,:) = new_prod.x_rand;
        
        x_rand2_l(ind, :) = new_prod.x_rand2;
        
        craft_l(ind) = new_prod.craft;
        
        k = k+1;
    end
    
end

if any(isnan(meanutility_l(:)))
    error('nan in demand');
end

if any(isnan(mc_mean_l(:)))
    error('nan in mc');
end


% the set of products that are fixed. data are already organized at the product level
if ~isempty(fixed_prod)
    w0 = [w0; fixed_prod.wp];
    p0 = [p0; fixed_prod.p0];
    
    carrier_ind_l = [carrier_ind_l; fixed_prod.retail_ind];
    oem_ind_l = [oem_ind_l; fixed_prod.owner_ind];
    
    if n_draws==1
        meanutility_l = [meanutility_l; fixed_prod.demand_resid];
        mc_mean_l = [mc_mean_l; fixed_prod.mc];
    else
        meanutility_l = [meanutility_l; fixed_prod.demand_resid+zeros(1, n_draws)];
        mc_mean_l = [mc_mean_l; fixed_prod.mc+zeros(1, n_draws)];
    end
    
    model_id_l = [model_id_l; fixed_prod.id];
    
    x_rand_l = [x_rand_l; fixed_prod.x_rand];
    
    x_rand2_l = [x_rand2_l; fixed_prod.x_rand2];
    
    craft_l = [craft_l; fixed_prod.craft];
    
    ind_endo_l = [ind_endo_l; true(size(fixed_prod.ind_endogenous))];
    
end

n_prod_l = size(x_rand_l, 1);

% compute profits
if n_prod_l>0
    oem_ind_l0 = oem_ind_l;

    p_sim = nan(n_prod_l, n_draws);
    w_sim = nan(n_prod_l, n_draws);
    share_sim = nan(n_prod_l, n_draws);
    carrier_profit_sim = nan(n_draws, 1);
    oem_profit_sim = nan(n_oem, n_draws);
    cons_surplus_sim = nan(n_draws, 1);
    cons_surplus_no_error_sim = nan(n_draws, 1);
    expmu_sim = cell(n_draws, 1);
    
    prod_profit_vec_ns = nan(size(p0, 1), n_draws);
    
    elas_ns = nan(size(p0, 1), size(p0, 1), n_draws);

    same_firm = false(n_prod_l, n_prod_l);
    for nc = 1 : size(oem_ind_l, 2)
        ind_nc = find(oem_ind_l(:, nc));
        if any(ind_nc)            
            same_firm(ind_nc, ind_nc) = true;
        end
        
    end
    
    for ns = 1 : n_draws
        [p_sim(:, ns), w_sim(:, ns), share_sim(:, ns), ...
            carrier_profit_sim(ns), oem_profit_sim(:, ns), ...
            cons_surplus_sim(ns), cons_surplus_no_error_sim(ns), ...
            expmu_sim{ns}, prod_profit_vec_ns(:, ns), elas_ns(:, :, ns)]...
            = equi_pwq(para_nonlinear, randomdraw, ...
            p0, w0,...
            x_rand_l, meanutility_l(:, ns), ind_endo_l, carrier_ind_l>0, oem_ind_l>0, ...
            ones(size(p0))*sim_pcoeff, ones(size(p0))*rand_income, income_norm, mc_mean_l(:, ns), setup, ...
            x_rand2_l, same_firm);
    end

    oem_ind_l = oem_ind_l0;
    p = mean(p_sim, 2);
    w = mean(w_sim, 2);
    share = mean(share_sim, 2);
    carrier_profit = mean(carrier_profit_sim, 2);
    oem_profit = mean(oem_profit_sim, 2);
    cons_surplus = mean(cons_surplus_sim);
    cons_surplus_no_error = mean(cons_surplus_no_error_sim);
    expmu = expmu_sim;
    prod_profit_vec = mean(prod_profit_vec_ns, 2);
    
    prod_profit_binary = prod_profit_vec;
else
    oem_ind_l = [];
    p = [];
    w = [];
    share = [];
    carrier_profit = [];
    oem_profit = zeros(n_oem, 1);

    oem_profit_sim = [];

    cons_surplus = [];
    cons_surplus_no_error = [];
    expmu = [];
    

    p_sim = [];
    w_sim = [];
    share_sim = []; 
    
    prod_profit_binary = []; 
    
    elas_ns = []; 
    
    craft_l = [];
end

if any(isnan(carrier_profit)) || any(isnan(oem_profit)) || any(isnan(cons_surplus))
    save('ws_nan.mat');
    error('nan in profit or consumer surplus');
end